Conditional Mean: \(0\)
Conditional Variance : \(\sigma_t^2\)
Daily Price of SP500 ETF (SPY) from Jan 02 2000 to Dec 31 2014
## Warning in if (!header) rlabp <- FALSE: the condition has length > 1 and
## only the first element will be used
## Warning in if (header) {: the condition has length > 1 and only the first
## element will be used
## ixDay Date Weekday X SPY.Open SPY.High SPY.Low SPY.Close
## 1 1 1/3/2000 Monday 1 148.25 148.25 143.88 145.44
## 2 2 1/4/2000 Tuesday 2 143.53 144.06 139.64 139.75
## 3 3 1/5/2000 Wednesday 3 139.94 141.53 137.25 140.00
## 4 4 1/6/2000 Thursday 4 139.62 141.50 137.75 137.75
## 5 5 1/7/2000 Friday 5 140.31 145.75 140.06 145.75
## 6 6 1/10/2000 Monday 6 146.25 146.91 145.03 146.25
## SPY.Volume SPY.Adjusted
## 1 8164300 110.33
## 2 8089800 106.01
## 3 12177900 106.20
## 4 6227200 104.50
## 5 8066500 110.57
## 6 5741700 110.94
layout(1)
# Fit GARCH model
library(fGarch)
Fit1 <- garchFit(~ garch(1,1), data=X1, cond.dist="norm", include.mean = FALSE, trace = FALSE)
## Fit1@fit$par #-- estimated parameters
## Fit1@residuals #-- this is not the garch residuals!!!! same as X1
## Fit1@sigma.t #-- estimated sig_t
## Fit1@residuals/Fit1@sigma.t #-- this is the (standardized) GARCH residuals
print(Fit1@fit$ics) #-- AIC and BIC are here## AIC BIC SIC HQIC
## -6.470319 -6.465360 -6.470321 -6.468556
plot(X1)
lines( 1.96*Fit1@sigma.t, type="l", col="red")
lines(-1.96*Fit1@sigma.t, type="l", col="red")## B-L test H0: the sereis is uncorrelated
## M-L test H0: the square of the sereis is uncorrelated
## J-B test H0: the sereis came from Normal distribution
## SD : Standard Deviation of the series
## BL15 BL20 BL25 ML15 ML20 JB SD
## [1,] 0 0 0 0.719 0.855 0 1.001
Starting with the definition of GARCH(1,1),
\[ \begin{align} \sigma_t^2 &= \omega + \alpha Y_{t-1}^2 + \beta \sigma^2_{t-1} \\\\ &= \omega + \alpha Y_{t-1}^2 + \beta \Big(\omega + \alpha Y_{t-2}^2 + \beta \sigma^2_{t-2} \Big)\\ \\ &= \omega + \beta \omega + \alpha Y_{t-1}^2 + \beta \alpha Y_{t-2}^2 + \beta^2 \sigma^2_{t-2} \\ \\ &= \omega + \beta \omega + \alpha Y_{t-1}^2 + \beta \alpha Y_{t-2}^2 + \beta^2 \Big(\omega + \alpha Y_{t-3}^2 + \beta \sigma^2_{t-3} \Big) \\ \\ &= \omega + \beta \omega + \beta^2 \omega + \alpha Y_{t-1}^2 + \beta \alpha Y_{t-2}^2 + \beta^2 \alpha Y_{t-3}^2 + \beta^3 \sigma^2_{t-3} \end{align} \]
Continuing, we get \[ \begin{align} &= \omega (1 + \beta + \beta^2 + \cdots ) + \alpha \sum_{i=0}^k \beta^i Y_{t-1-i}^2 + \beta^{k+1} \sigma^2_{t-1-k} \\\\ &= \frac{\omega}{1-\beta} + \alpha \sum_{i=0}^\infty \beta^i Y_{t-1-i}^2 \end{align} \]
That means if we keep going, we can write \[ \sigma_t^2 \hspace3mm = \hspace3mm \frac{\omega}{1-\beta} + \alpha \sum_{i=0}^\infty \beta^i Y_{t-1-i}^2 \]
We will use the truncated and estimated version of this, \[ \hat \sigma_t^2 \hspace3mm = \hspace3mm \frac{\hat \omega}{1- \hat \beta} \hspace3mm + \hspace3mm \hat \alpha \sum_{i=0}^{t-1} \hat \beta^i Y_{t-1-i}^2 \]
GARCH(1,1) says, \[ Y_t \hspace3mm = \hspace3mm \sigma_t e_t \]
Using observation \(\{Y_1, \ldots, Y_n\}\), the residuals are \[ \hat e_t = Y_t/ \hat \sigma_t \]
\[ \begin{align} Y_t &= \sigma_t e_t \hspace{10mm} e_t \sim_{iid} N(0,1)\\\\ \sigma_t^2 &= \omega + \alpha Y_{t-1}^2 + \beta \sigma^2_{t-1} \end{align} \]
[Conditonal distribution] \(=\) [Distribution of \(e_t\)]
Guessing correct distribution for \(e_t\) is very important in GARCH parameter estimation.
x <- seq(-10,10, .1)
#- Standardized t-distribution
plot(x, dstd(x, mean = 0, sd = 1, nu = 5), type="l" )
lines(x, dnorm(x, mean = 0, sd = 1), col="red" ) #- Generalized Error Distribution
plot(x, dged(x, mean = 0, sd = 1, nu = 4), type="l" )
lines(x, dnorm(x, mean = 0, sd=1) , col="red" ) #- Skewed-Standardized t-distribution
plot( x, dsstd(x, mean = 0, sd = 1, nu = 5, xi = 1.5), type="l" )
lines(x, dnorm(x, mean = 0, sd = 1), col="red" ) #- Skewed-Generalized Error Distribution
plot(x, dsged(x, mean = 0, sd = 1, nu = 5, xi = 1.5) , type="l" )
lines(x, dnorm(x, mean = 0, sd = 1), col="red" )True parameter (.024, .1, .8)
True conditional distribution: Normal Estimated using: Normal
True conditional distribution: std(5) Estimated using: Normal
True conditional distribution: sged(skew=.7, shape=1.45) Estimated using: Normal
#--- Simulation showing how error distribution affect parameter estimation in GARCH(1,1) ---
n <- 1000
itt <- 1000
theta <- c(.024, .1, .8)
Param11 <- matrix(rep(0,3*itt), itt, 3)
for (i in 1:itt) {
spec <- garchSpec(model = list(omega=theta[1], alpha=theta[2], beta=theta[3]), cond.dist="norm")
x <- garchSim(spec, n = n, extended=FALSE)
est2 <- garchFit(~ garch(1,1), data=x, cond.dist="norm", include.mean = FALSE, trace = FALSE)
Param11[i,] <- coef(est2)
}
Param12 <- matrix(rep(0,3*itt), itt, 3)
for (i in 1:itt) {
spec <- garchSpec(model = list(shape=5, omega=theta[1], alpha=theta[2], beta=theta[3]), cond.dist="std")
x <- garchSim(spec, n = n, extended=FALSE)
est2 <- garchFit(~ garch(1,1), data=x, cond.dist="norm", include.mean = FALSE, trace = FALSE)
Param12[i,] <- coef(est2)
}
Param13 <- matrix(rep(0,3*itt), itt, 3)
for (i in 1:itt) {
spec <- garchSpec(model = list(skew=.9, shape=1.45, omega=theta[1], alpha=theta[2], beta=theta[3]), cond.dist="sged")
x <- garchSim(spec, n = n, extended=FALSE)
est2 <- garchFit(~ garch(1,1), data=x, cond.dist="norm", include.mean = FALSE, trace = FALSE)
Param13[i,] <- coef(est2)
}
layout(matrix(1:9, 3, 3, byrow=T))
hist(Param11[,1], 10, xlim=c(0,.5)); hist(Param11[,2], 10, xlim=c(0,.5)); hist(Param11[,3], xlim=c(0.4,1), breaks=10)
hist(Param12[,1], 10, xlim=c(0,.5)); hist(Param12[,2], 10, xlim=c(0,.5)); hist(Param12[,3], xlim=c(0.4,1), breaks=10)
hist(Param13[,1], 10, xlim=c(0,.5)); hist(Param13[,2], 10, xlim=c(0,.5)); hist(Param13[,3], xlim=c(0.4,1), breaks=10)
P11 <- (Param11 - t(matrix(theta, 3, itt) ))^2
P12 <- (Param12 - t(matrix(theta, 3, itt) ))^2
P13 <- (Param13 - t(matrix(theta, 3, itt) ))^2
apply(P11, 2, mean)
apply(P12, 2, mean)
apply(P13, 2, mean)
apply(P11, 2, mean)/apply(P11, 2, mean)
apply(P12, 2, mean)/apply(P11, 2, mean)
apply(P13, 2, mean)/apply(P11, 2, mean)MSE1/MSE1 #- Relative MSE
[1] 1 1 1
MSE2/MSE1
[1] 1.584734 3.396001 1.831015
MSE3/MSE1
[1] 1.198293 1.613116 1.355613GARCH(1,1) is weakly stationary if \[ E( \log(\beta + \alpha e_t^2)) < 0 \]
This is satisfied if \[ \alpha + \beta < 1. \]
library(quantmod)
library(fGarch)
source("http://gozips.uakron.edu/~nmimoto/477/TS_R-90.txt")
theta = c(.5, .7, .2)
spec1 <- garchSpec(model = list(omega=theta[1], alpha=theta[2], beta=theta[3]), cond.dist="norm")
#- Generate GARCH
set.seed(66343)
X <- garchSim(spec1, n = 1000, extended=T) #extended=T: now Y has 3 columns: garch sigma eps## Warning in if (extended) ans else ans[, "garch"]: the condition has length
## > 1 and only the first element will be used
X <- as.xts(X)
X <- xts(X, order.by=as.Date(index(X)))
Y <- X[,1] #- simulated GARCH
sig.true <- X[,2] #- this is true sigma
plot(Y, main="Simulated GARCH") #- Estimate GARCH
out1 <- garchFit(~ garch(1,1), data=Y, cond.dist="norm", include.mean = FALSE, trace = FALSE)
sig.estim <- xts(out1@sigma.t, order.by=index(Y))
#- compare true sigma vs estimated sigma
head( cbind(sig.true, sig.estim) )## sigma ..2
## 2015-07-14 0.9844680 1.7773189
## 2015-07-15 0.8719402 1.0198346
## 2015-07-16 0.8528739 0.9322072
## 2015-07-17 0.9465185 1.0167695
## 2015-07-18 0.9430606 1.0033031
## 2015-07-19 1.1075998 1.1624516
#- Compute the sigma for tomorrow (forecast sigma_t)
w <- out1@fit$par[1]
a <- out1@fit$par[2]
b <- out1@fit$par[3]
sig.pred <- xts( sqrt( w + a*Y[1000]^2 + b*out1@sigma.t[1000]^2 ), order.by=index(Y[1000])+1)
#- This does the same thing (10-day forecast sigma_t)
sig.pred2 <- predict(out1)
sig.pred2 ## meanForecast meanError standardDeviation
## 1 0 1.061806 1.061806
## 2 0 1.267963 1.267963
## 3 0 1.411732 1.411732
## 4 0 1.517237 1.517237
## 5 0 1.596835 1.596835
## 6 0 1.657922 1.657922
## 7 0 1.705340 1.705340
## 8 0 1.742442 1.742442
## 9 0 1.771642 1.771642
## 10 0 1.794720 1.794720
## [1] 1.061806
#- Rolling 1-day prediction of sig.t
sig.pred1 <- numeric(0)
date.stp <- numeric(0)
for (i in 1:750){
T <- Y[(1:250)+(i-1)]
n <- length(T)
out1 <- garchFit(~ garch(1,1), data=T, cond.dist="norm", include.mean = FALSE, trace = FALSE)
sig.pred1[i] <- predict(out1)[1,3] #- sig.t prediction
date.stp[i] <- index(T[n])+1
}
sig.pred <- xts(sig.pred1, order.by=as.Date(date.stp))
#- Plot and compare estimated vs rolling predicted
plot(cbind(sig.pred, sig.estim, sig.true),
col=c("blue", "red", "black"),
main="sigma.t: True (Blk) vs Estim (Red) vs Pred(Blue)") plot(cbind(sig.pred, sig.estim, sig.true)['2016'],
col=c("blue", "red", "black"),
main="sigma.t: True (Blk) vs Estim (Red) vs Pred(Blue)") plot(cbind(sig.pred, sig.estim, sig.true)['2016-04'],
col=c("blue", "red", "black"),
main="sigma.t: True (Blk) vs Estim (Red) vs Pred(Blue)") ts.plot(cbind(Y, 1.96*sig.true, -1.96*sig.true, 1.96*sig.pred, -1.96*sig.pred),
col=c("black", "black", "black", "blue", "blue")) ts.plot(cbind(Y, 1.96*sig.true, -1.96*sig.true, 1.96*sig.pred, -1.96*sig.pred)['2016-06::2016-12'],
col=c("black", "black", "black", "blue", "blue"))## [1] "GSPC"
#- we model this with GARCH (this time we don't know the true sigma.t)
Y <- diff(log(SPX["2013::2016"]))[-1]
length(Y)## [1] 1007
#- Estimate GARCH
out1 <- garchFit(~ garch(1,1), data=Y, cond.dist="norm", include.mean = FALSE, trace = FALSE)
sig.estim <- xts(out1@sigma.t, order.by=index(Y))
#- Rolling 1-day prediction of sig.t
sig.pred1 <- numeric(0)
date.stp <- numeric(0)
for (i in 1:750){
T <- Y[(1:250)+(i-1)]
out1 <- garchFit(~ garch(1,1), data=T, cond.dist="norm", include.mean = FALSE, trace = FALSE)
sig.pred1[i] <- predict(out1)[1,3] #- sig.t prediction
date.stp[i] <- index(T[length(T)]+1)
}
sig.pred <- xts(sig.pred1, order.by=as.Date(date.stp))
#- Compare estimated sigma vs predicted sigma -
plot(cbind(sig.pred, sig.estim), col=c("blue", "red")) #- Compare daily CI from two sigmas
plot(as.numeric(Y["2016"]), type="h")
lines(as.numeric( 1.65*sig.estim["2016"]), col="red")
lines(as.numeric(-1.65*sig.estim["2016"]), col="red")
lines(as.numeric( 1.65*sig.pred["2016"]), col="blue")
lines(as.numeric(-1.65*sig.pred["2016"]), col="blue")Cryer : CREF stock values
Shumway : NYSE us GNP
Cowpertwait : SP500, Southern hem Temp.